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^ '. Abstract 



We study the effect of confineinent in the dynamical behavior of a core-softened fluid. The fluid 
is modeled as a two length scales potential. This potential in the bulk reproduces the anomalous 
behavior observed in the density and in the diffusion of liquid water. A series of NpT Molecular 
Dynamics simulations for this two length scales fluid confined in a nanotube were performed. We 
obtain that the diffusion coefficient increases with the increase of the nanotube radius for wide 
channels as expected for normal fluids. However, for narrow channels, the confinement shows 
an enhancement in the diffusion coefficient when the nanotube radius decreases. This behavior, 
observed for water, is explained in the framework of the two length scales potential. 



I. INTRODUCTION 

The dynamic behavior of fluids in the bulk is characterized by transport properties such 
as the diffusion coefficient. In simple liquids they are governed by the molecular interactions 
which can be effectively treated as pair potentials. In complex fluids, such as water, the 
dynamics of the systems is governed by network-forming interactions. The potentials in- 
volved are in general highly anisotropic and the transport properties tend to exhibit unusual 
behaviors. 

Under confinement even normal liquids have an unusual behavior, very different from 
the physical properties observed in bulk. The competition between surface effects and the 
confinement can induce a dramatic change in the transport properties of the fluid inside the 
channel^"-. 

Bulk water is anomalous in many of its characteristics. The maximum in water's density 
is a well-known anomaly but there are many others. The self-diffusion coefficient at fixed 
temperature for a normal liquid decreases under compression, while in liquid water it in- 
creases with the increase of pressure. In bulk water this is due to the hydrogen bonds that 
are created and destroyed making particles to move from one neighbor to another neighbor. 

Notwithstanding its molecular simplicity water is quite hard to be modeled. The reason 
behind this difficulty is the presence of the hydrogen bonds, a non symmetric charge dis- 
tribution and polarizability of the molecule that are density and temperature dependents. 
Consequently, there are more than twenty-five (bulk) water models for computational simu- 
lation - empirical potentials - in which each of them give a different dipole moment, dielectric 
and self-diffusion constants, average configurational energy, density maximum, and expan- 
sion coefficient. More specifically, the maximum of density is found experimentally to be 
at T = 4 °C (for pressure P = 1 atm) while such models give values ranging from -45 °C 
(SPC model) up to 25 °C (P0L5/TZ model). The TIP5P water model was built to match 
the 4 °C experimental result, but it fails in many other aspects.- Despite these limitations, 
these models have been used to understand the transport properties and phase transitions 
of confined water-"—. The results give a qualitative comparison with experiments without 
providing a complete understanding of the origin of the anomalies.— The majority of the 
molecular water models are conceived focusing on accurately describe the hydrogen bonds 
and charge distributions since many of the water uncommon properties are believed to come 



from its highly directional interactions. Examples are solvation and properties which depend 
on polarization. On the other hand, the literature have many examples in which systems 
with absence of anisotropic interactions still may present some of the water features^"—. 
Some of its anomalous behavior may come from purely volumetric effects, which particularly 
is our focus in this work. 

For confined systems, where water molecules interact in nanoscale distances, first princi- 
ples simulations would be the appropriated tool for numerical comparison with experimental 
data. This procedure, however, has limitations. Even for confined water systems, in which 
the sizes involved are much smaller than that ones found in bulk cases, thousands of atoms 
are necessary for attacking typical problems along with millions of simulation steps. In this 
sense, turns out that in the majority of cases ab initio techniques become impracticable for 
dealing with such computational demanding systems. 

Given the limitations of the full water models and the computational costs of the ab 
initio simulations, classical effective empirical potentials became the simplest framework to 
understand the physics behind the anomalies of bulk water. From the desire of constructing 
a simple two-body potential capable of describing the anomalous behavior of bulk water, 
a number of models have been developed^"—. Despite their simplicity, such models had 
successfully reproduced the thermodynamic, dynamic, and structural anomalous behavior 
present in bulk liquid water. They also predict the existence of a second critical point 
hypothesized by Poole and collaborators^^. This suggests that some of the unusual properties 
observed in water can be quite universal and possibly present in other systems. 

In the case of confined water a number of attempts have been made to understand its 
thermodynamic and dynamic properties. For the confinement media, nanotubes have been 
widely used for mimicking water confined into live organisms and as building blocks for 
technological applications, as desalination of water'^°. Also, nanotubes can be used for drug 
delivery since they resemble biological ion channels^. In addition, confinement in nanopores 
and nanotubes have been used also to avoid spontaneous water crystallization below the 
melting point in an attempt to observe its hypothetical second critical point^"— . Simulations 
employing some of the discussed molecular models for water, namely SPC/E, TIP4P-EW 
and ST2, confined in nanoscale channels exhibit two complementary effects: the melting 
temperature of the fluid at the center of the channel decreases and water crystallizes at the 
channel surface^"—. In addition to these thermodynamic properties, the mobility properties 



of confined water also exhibit an unusual behavior. Experiments show an enhancement of 
the mobility orders of magnitude higher than what is predicted by the flow theories^i^. 
Simulations^"— show an increase in the enhancement rate below a certain threshold radius. 
Similarly the self-diffusion coefficient, D, obtained through molecular dynamic simulations 
for atomistic models, below a certain radius increases with with decreasing radius^^i^"— . 

Besides the thermodynamic and dynamic unusual properties of confined water, the struc- 
ture also presents an interesting behavior. The water structure inside larger nanotubes 
exhibits a layered structure, while in narrow nanotubes a single file is observed. The layered 
water molecules can be found in a spiral-like chain,— in a hexagonal structure for (6,6) car- 
bon nanotubes (CNT),— a octagonal water-shell structure for a (9,9) CNT or a octagonal 
water-shell structure with a central water chain for a (10,10) CNT,— i^ and others different 
structures^. 

The presence of layering effects is also controversial. Molecular Dynamics (MD) simula- 
tions results by Wang et alM do not show any obvious ordered water structure for (9,9) or 
(10,10) CNTs, unlike the works of Koles et aP^^. This difference in findings can be caused 
by the infiuence of the water model used and corresponding Lennard- Jones parameters^. 
The molecular structure of water confined in nanotubes and the diffusion can be very dif- 
ferent depending on the chosen water model to perform the simulations^^. This difference 
arises as a consequence of the fact that the water models used in classical all-atoms MD 
simulations, like SPC/E, TIP3P, TIP4P, etc, are parametrized for bulk simulations, and for 
reproduce only few aspects of real water. So this models may show errors to describe the 
correct water behavior under strongly confinement, like inside nanotubes. 

The behavior of the diffusion coefficient with the nanotube radius is still under debate. 
Some simulations show a monotonical decrease of D with decreasing nanotube radius^"— 
while other simulations indicate the presence of a minimunt^^i^"— . It has been suggested 
that the length of the nanotube and the length of the simulation would be responsible for 
the different results^'^'^ and that friction should play also a relevant role^. 

Our core-softened model introduced to study bulk system does not have any directionality 
and therefore it is not water. However, it does exhibit the density, the diffusion and the 
response functions anomalies observed in water. This suggests that some of the anomalous 
properties that are attributed directionality of water can be found in spherical symmetry 
systems. Likewise, here we propose that also the minimum in the diffusion coefficient as the 



nanotube radius is decreased observed for water can be also found in spherical symmetric 
systems. 

In order to check our hypothesis, we model a water-like fluid into a nanotube using a 
core-softened potential. We test if this model is capable to capture the increase in the 
diffusion coefficient when the channel radius decreases. Next, we verify if the layering and 
the structure formed inside the channel, observed in some classical models for confined liquid 
water, has an universal feature or if it is just a consequence of the specific confining surface 
and model details. 

The paper is organized as follows. The nanotube, water-like fiuid model and the sim- 
ulational details are presented in Sec. [ITl Our results are discussed in Sec. IIIIl and the 
conclusions and summary are presented in Sec. IIVI 

II. THE MODEL AND THE SIMULATION DETAILS 
A. The Model 




FIG. 1: Interaction potential between water-like particles. 



The water-like fiuid is modeled as point particles with effective diameter a and mass m, 
interacting through the three dimensional core-softened potentiay^*^ 
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The first term in Eq. ([T]) is the standard Lennard- Jones (LJ) 12-6 potential^ and the second 
term is a Gaussian centered in ro/a, with depth u^e and width ca. For uq = 5.0, c = 1.0 and 
^o/c" = 0.7 this equation represents a two length scale potential, with one scale at Tjj ^ 1.2cr, 
when the force has a local minimum, and the other scale at Tjj ^ 2a, where the fraction 
of imaginary modes has a local minimum^, de Oliveira et al.—^ obtained the pressure- 
temperature phase diagram of this system and showed that it exhibits thermodynamic, 
dynamic and structural anomalies similar to the anomalies present in water— i^. 

Here we study the dynamic behavior of this water-like model confined in a nanotube 
connected to two reservoirs. The nanotube-reservoir setup is illustrated in Fig. [2j The 
simulation box is a parallelepiped with dimensions L^. x Ly x L^. 

Two fluctuating walls, A in left and B in right, are placed in the limits of the x-direction of 
the simulation box. The walls are allowed to move in order to maintain the pressure constant 
in the reservoirs. The sizes Ly and of L^ depend on the effective nanotube radius, a, and 
they are defined hj Ly = Lz = L = 2a + 6a. The initial size L^ is given by Lx = QLc, where 
Lc is the tube length. The nanotube structure was constructed as a wrapped hexagonal 
lattice sheet of point particles whose diameter is o"nt = ct- The nanotube interacts with the 
water-like particles through the Weeks- Chandler- Andersen (WCA) potential^° given by 
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where Ui^^ir) is the standard LJ potential. The cutoff distance for this interaction is Tc = 
2}l^aij, where aij = (ai + aj)/2 is the center-to-center distance between the fluid particle i 
and the nanotube particle j. 

B. The simulation details 

The properties of the system were evaluated with simulations at constant number of 
particles, pressure and temperature {NpT ensemble). The Andersen thermostat^ with 
collision frequency uSt = 0.01 was used to maintain the temperature flxed. The pressure in 
both reservoirs was flxed using the Lupowski and van Smol method of fluctuating conflning 
walls^^'^^. These fluctuating walls act like pistons in the system where a constant force 
controls the pressure in the x-direction. This lead us to rewrite the resulting force in a 
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FIG. 2: Schematic depiction of the simulation box with the nanotube, reservoirs and fluctuating 
walls. The cylindrical channel in the center has radius a and length Lc- The reservoirs have height 
L. 

water-like particle as 

Fr = -VU^j + Fi^A{nA) + Fi^B{riB) , (3) 

— * 

where Fi^j indicates the interaction between the particle i and the piston j. These forces were 
calculated from a WCA potential similar to Eq. (J2]), however considering the x-projection 
of the distance between one particle in the bulk and the piston position. 
The equation of motion for the pistons are 
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m^aB = pS^HB - ^ Fi^Biris) , (5) 

where m^ is the piston mass, p the desired pressure in the system, 5*^ is the piston area and 
Ha is a unitary vector in positive x-direction, while hb is a negative unitary vector. Both 
pistons (A and B) have mass m^ = m = 1, width o"^ = a and area equal to 5*^ = L^. The 
Andersen thermostat is also applied to the pistons to ensure the temperature control. The 
values of pressure and temperature were chosen avoiding the density anomaly and the solid 
state regionai^ii^. 

For simplicity, we assume that the nanotube atoms are fixed (i.e., not time integrated) 
during the simulation. The reduced quantities are defined as usual. 
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for the channel radius, density of particles, time and temperature, respectively, and 

p* = - — and D* = ^ ' ' — (7) 

e a 

for the pressure and diffusion coefficient, respectively. Periodic boundary conditions were 
applied in the y and z directions. The equations of motion for the particles of the fluid were 
integrated using the velocity Verlet algorithm, with a time step (5t* = 0.005. The fluid-fluid 
interaction, Eq. ([1]), has a cutoff radius r*^^ = 3.5. The nanotube radius was varied from 
a* = 1.25 to a* = 10.0, and the number of fluid particles in the simulations varies from 
500 to 3500. The number of particles were chosen considering that the nanotube would be 
filled with the fiuid and that we would have in the reservoirs the same properties evaluated 
in previous NVT simulation for the non-confined casei^^^^. For all values of radius the 
nanotube length was defined as L* = 20. 

Five independent runs were performed to evaluate the properties of the fiuid inside the 
nanotube. For each simulation run half of fiuid particles was initially placed into each 
reservoir. We performed 5 x 10^ steps to equilibrate the system followed by 5 x 10® steps for 
the results production stage. The equilibration time was taken in order to ensure that the 
nanotube became filled with water-like particles as well as the pistons reached the equilibrium 
position for a given pressure. 

For calculating the axial diffusion coefficient, D^, we computed the axial mean square 
displacement (MSD) namely 

{[x{t) - x{to)]') = {Ax{tf) = 2Dt\ (8) 

where x{to) and x{t) denote the axial coordinate of the confined water-like molecule at a 
time to and a later time t, respectively. The diffusion coefficient Dx is then obtained from 

t— >oo Zt 

Depending on the scaling law between Ax^ and t in the limit t ^^ oo, different diffusion 
mechanisms can be identified: a = 0.5 identifies a single file regime^, a = 1.0 stands for a 
Fickian diffusion whereas a = 2.0 refers to a ballistic diffusion^i'^'^i^i^. 

III. RESULTS AND DISCUSSION 

First, we checked which is the diffusive regime of our system for different channel ra- 
dius. Fig. |3] illustrates the axial mean square displacement versus time for channel radius 

8 



a* = 1.25, 1.5, 2.0, 5.0, and 10.0 at T* = 0.25 and p* = 0.7. For simple LJ confined fiuids 
for very narrow channels a single-file diffusion regime is found {a = 0.5)^. Instead, in 
our model Fickian diffusion was observed for all channel radius, {a = 1.0). This result is 
in agreement with the diffusion coefficient at the center of nanotubes observed for water 
models SPC/EP^^s , SPG^, TIP3P^ and TIP4P^. In the case of water the Fickian dif- 
fusion has been interpreted as a combination of highly coordinated like in a ballistic mode 
and a mono dimensional configuration like in a single-line diffusiort^'^. The presence of a 
highly coordinated structure was observed in water in very narrow nanotubes^ and between 
plates^, confirming the interpretation of Mukherjee et al.— . In our case, the particles also 
combine a strongly correlated motion with a single-line diffusion. As we are going to see in 
detail below, the molecules are arranged in shells. The motion inside each cell is correlated 
and slow while the motion between shells is faster. The compromise between these two 
mechanisms lead to a Fickian diffusion. 

It is important to stress that the major difference between a simple LJ fluid and our 
core-softened model is the presence of two-scales in the potential (absent in the LJ fluid). 
This enforces our conjecture that two scales in the interatomic potential plays an important 
role in the appearance of water-like features. 
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FIG. 3: Axial mean square displacement versus time for channel radius a* = 1.25, 1.5, 2.0, 5.0, 
and 10.0. 



Next, we tested if the diffusion through the channel obeys the mean-fleld-like Knudsen 
equation, i.e., if the diffusion coefficient is proportional to the channel radius. Fig. H] il- 
lustrates the diffusion coefficient, D, versus channel radius, a*, for fixed T* = 0.25 and 
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FIG. 4: Diffusion coefficient inside nanotube, D^^nt, in units of non-confined diffusion, Dx^buik, for 
different nanotube radius. Tfie error bars are smaller than the data point. The dotted line is a 
guide to the eye. 

p* = 0.7. We see from this figure that a critical channel radius a* exists where the derivative 
of the D{a) curve is zero. For a* > a^ = 2.0 the diffusion coefficient presents the expected 
behavior of increasing with a*. For large channel radius the growth is linear as predicted 
by the Knudsen equation. For a* < 2.0, on the other hand, we observed that D decreases 
with increasing a*, which can not be explained by the Knudsen mean- field approach. At 
a* = al = 2.0 particles are virtually immobilized, i.e., D ^ 0. 

Studies for SPC/E^^'^, TIP4P-EW^^ and SPC^^'^e g^^^^ ^j^^t the diffusion coefficient 
increases with the channel radius, a. Simulations for SPC/E^i^ and TIP4P-EW— 1^ also 
show the decrease of the diffusion coefficient with the increase of the channel radius for 
a < ttc- This anomalous region is captured by our model and it is not observed in our 
simulations for LJ confined fluids (details not shown here for simplicity). 

For SPC/E and TIP4P-EW potentials used for confined water, the number of neighbors 
and the number of hydrogen atoms differ from those numbers in the bulk phase. In these 
models the different slope in the D{a) function, i.e., positive for a > ac and negative for a < 
Qci are attributed to a competition between two effects: the confinement and the nanoscale 
surface. For the a > ac case, D decreases for decreasing a because of the confinement. This 
is not hard to understand since decreasing a allows less space for particles to move^i^"— . 
Increasing confinement leads to surface effects becoming more important. In the water case 
hydrogen bonds from the surface are depleted and molecules become more mobile. This 



10 



would explain why D increases for decreasing a below a,Mi^^. The behavior of D in water, 
therefore, can be explained by minimizing the free energy. Water molecules have a gain in 
rotational entropy when they are come inside a narrow carbon channel^>^. This gain gain 
in entropy compensates the lost in elthalpy due to the reduction of number of hydrogen 
bondsi^"^. 

Interestingly though our system does not have hydrogen bonds, thus it is not subject 
to the competition between hydrogen bonds depletion and diminishing space available for 
particles to move. Therefore, what would be the mechanism behind our non-monotonic 
curve D{a) shown in Fig. HP 
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FIG. 5: Radial density profile for different values of radius: (A) a* 
a* = 2.0, (D) a* = 4.0, (E) a* = 7.0 and (F) a* = 10.0. 



1.25, (B) a* = 1.5 (C) 



The behavior of the diffusion coefficient in our model can be understood by two com- 
plementary ways. First, by examining the density profile inside the nanotube. The density 
distribution is computed in cylindrical coordinates, r"^ = y"^ + z"^, where r = is the center 
of the channel. Fig. [5] illustrates the radial density profile versus r* for the channel radii 
a* = 1.25, 1.5, 2.0, 4.0, 7.0, and 10.0. In all analyzed cases layering is observed. Axial layers 
are also observed in simulations for the SPC/E and TIP4P-EW models for water confined 
in nanotubes^'^'^'^'^. In the last cases, the presence of layering is attributed to the hy- 
drogen bonds and surface effects. In our model the presence of layering comes as a result 
of the competition between particle-particle and particle-wall interactions. The potential 
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FIG. 6: Snapshots of the system for (A) a* = 2.0 and (B) a* = 1.5. 

illustrated in Fig. [1] favors particles to be at least at r*p = 2.0 apart, while the hydrophobic 
walls push particles away to a distance of at least r*^ = 2^/^. Consequently, for a* = 1.25, 
1.5, 2.0, 4.0, and 7.0 a number of layers equal to 1, 1, 2, 4 and 6 are formed. For a* > 10 a 
continuous distribution emerge. For a* > 2.0 the system form layers arranged in distances 
that minimize the potential energy. This indicates that for larger diameters the enthalpic 
contribution for the free energy dominates over the entropic contribution. For a* < 2.0 
there is only one layer of particles. Since the particle-wall interaction is purely repulsive, 
particles advance over the wall repulsive region moving "free" in the radial direction and 
entropy increases and the mobility rises. For a* < 2.0 the entropic contribution dominates 
over the enthalpic contribution. 

The decrease in the diffusion coefficient as the channel radius is decreased for a* > a* = 
2.0 is associated with the layers formation and particularly with the correlation between 
particles in different layers that try to move without changing the layer to layer distance. As 
the number of layers increase for a* > a* = 2.0 fluctuations allow particles to move faster. 
At a* = a* = 2.0 the diffusion reaches a minimum and the system assumes a crystal-like 
configuration as illustrated in Figured Therefore, the confinement leads the fiuid to a solid- 
like state, even at values of temperature and pressure far from solid state phase^i^i^i^i^. 

In order to check what happens as the channel radius is decreased further, snapshots of 
the system for a* = 2.0 and a* = 1.5 are shown in Fig. [61 For a* = 2.0 particles form 
two layers while for a* = 1.5 a single layer is observed. The correlations that immobilizes 
particles at the two layers structure disappear as the single layer is formed and particles can 
diffuse faster for a* < 2.0 by moving from the close packing (at r* = 1.0) to the minimum 
of the fraction of imaginary modes (at r* = 2.0). 

Besides the layering in the radial direction, particles also change their structure in the 
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FIG. 7: Total density inside the nanotube as function of the tube radius for T* = 0.25 and p* = 0.7 
. The error bars are smaher than the data point. The hne is a guide to the eye. 



axial direction. Fig. [7] shows the total density inside the nanotube as function of the radius. 
For a* < 4.0 the density of the confined system increases with the decrease of a*. This 
result is qualitatively the same observed for recent simulations of SPC/E model of water 
in a nanotube-reservoirs system^. This can be explained by a change in the axial distance 
between particles. Fig. [8] for a* = 1.5 illustrates that for small radius the preferential axial 
distance is the shoulder scale, x* = 1.0, while in the bulk and for larger radius the preferential 
axial distance is x* ~ 2.0. 
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IV. CONCLUSION 

This paper explores the connection between the surface interaction, confinement and the 
presence of two length scales in the diffusion of a fluid in narrow channels. The water-like 
fluid was modeled using a spherically symmetric two length potential, and the conflning 
channel is modeled as hard spheres. Our system shows an enhancement of the diffusion co- 
efficient with the decrease of the channel radius for a channel radius below a certain critical 
value. This effect arises from the competition between the confinement, that accommodate 
particles at the lower energy length scale, and the surface interaction, that pushes parti- 
cles away from the surface generating correlated layers. For wider channels the layers are 
accomadated minimizing the potential energy forming an organized structure. For narrow 
channels the particles form a single layer that move more freely advancing over the wall 
repulsion and therefore increasing the diffusion. The mechanism is quite similar to the one 
proposed for water. 

In addition we found that below a certain channel radius the density inside the channel 
is larger than the bulk density. This is explained on basis of the length scales competition 
and it is in agreement with simulations. Our results indicate that the presence of minimum 
in the diffusion coefficient with the decrease of the channel radius is not a property solenely 
due to directional systems such as water but can be seen also in spherical symmetric systems 
or in system in which directionality would not play a relevant role. 
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